Referencia de la imagen: “The TibA Adhesin/Invasin from Enterotoxigenic Escherichia coli Is Self Recognizing and Induces Bacterial Aggregation and Biofilm Formation”
Asumimos que los puntos pueden tener cualquier localización en la ventana, en nuestro caso no hay ninguna limitación al respecto. El objetivo es estimar los parámetros de distribución de las diferentes bacterias.
En los datos tenemos dos tipos de poblaciones bacterianas (marks = 2).
#install.packages("imager")
#install.packages('spatstat')
library(imager)
library(spatstat)
library(cowplot)
library(ggplot2)
image <- load.image('./Bacterial_interaction.png')
plot(image, main = "Original image")
image <- resize(image,round(width(image)/4),round(height(image)/4))
plot(image,main="Low resolution")
image
## Image. Width: 112 pix Height: 112 pix Depth: 1 Colour channels: 4
Pasar la imagen a un dataFrame.
image.df <- as.data.frame(image)
dim(image.df)
## [1] 50176 4
Selección de los puntos con mayor intensidad de acuerdo al color verde y rojo (RGB). Con esta técnica algunos puntos de la imagen original pueden quedar excluidos. Se han seleccionado aquellos que en el correspondiente canal tengan más de un 50% de intensidad y el resto de canales menos del 0.1. Se ha realizado un filtrado de los valores que en el canal 3 (azul) tiene una intensidad superior a 0.3
image.df$quality <- image.df$cc
image.df[image.df$cc == 3 & image.df$value > 0.2, ]$quality = 0
spatial.points <- data.frame(x = image.df[image.df$value > 0.5 & (image.df$cc == 2 | image.df$cc == 1) & image.df$quality != 0, ]$x,
y = image.df[image.df$value > 0.5 & (image.df$cc == 2 | image.df$cc == 1) & image.df$quality != 0, ]$y,
strain = image.df[image.df$value > 0.5 & (image.df$cc == 2 | image.df$cc == 1) & image.df$quality != 0, ]$cc)
p <- ggplot(spatial.points, aes(x=x, y=y)) +
geom_point(aes(color = as.factor(strain))) + theme_minimal() + theme(legend.title = element_blank()) +
scale_color_manual(values = c("#FF0000", "#04FF00"))
p
Carga de los puntos con una ventana de dimensiones 112 x 112
pattern <- ppp(spatial.points$x, spatial.points$y, c(0, 112), c(0, 112))
plot(pattern, main = NULL)
summary(pattern)
## Planar point pattern: 716 points
## Average intensity 0.05707908 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are integers
## i.e. rounded to the nearest unit
##
## Window: rectangle = [0, 112] x [0, 112] units
## Window area = 12544 square units
pattern <- rescale(pattern, 10) # Rescale to lower distance
summary(pattern)
## Planar point pattern: 716 points
## Average intensity 5.707908 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
##
## Window: rectangle = [0, 11.2] x [0, 11.2] units
## Window area = 125.44 square units
intensity(pattern)
## [1] 5.707908
Error estandar
sqrt(intensity(pattern)/area(Window(pattern)))
## [1] 0.2133145
Q <- quadratcount(pattern, nx=6, ny=6)
Q
## x
## y [0,1.87) [1.87,3.73) [3.73,5.6) [5.6,7.47) [7.47,9.33)
## [9.33,11.2] 17 32 24 21 6
## [7.47,9.33) 21 29 15 18 12
## [5.6,7.47) 9 12 27 9 10
## [3.73,5.6) 28 13 19 21 31
## [1.87,3.73) 18 17 30 10 10
## [0,1.87) 21 27 21 21 25
## x
## y [9.33,11.2]
## [9.33,11.2] 33
## [7.47,9.33) 21
## [5.6,7.47) 17
## [3.73,5.6) 18
## [1.87,3.73) 24
## [0,1.87) 29
plot(intensity(Q, image = TRUE), main = "Intensity per quadrant GFP + RFP")
La hipótesis nula es que el patrón se distribuye de forma homogénea. En nuestro caso podemos rechazar dicha hipótesis de forma que aceptamos que nuestros datos no son homogéneos.
quadrat.test(pattern, nx=6, ny=6)
##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern
## X2 = 95.609, df = 35, p-value = 3.13e-07
## alternative hypothesis: two.sided
##
## Quadrats: 6 by 6 grid of tiles
Volviendo a realizar el test con un número menor de cuadrantes tenemos que se acepta la hipótesis nula.
quadrat.test(pattern, nx = 5, ny = 5)
##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern
## X2 = 47.617, df = 24, p-value = 0.005628
## alternative hypothesis: two.sided
##
## Quadrats: 5 by 5 grid of tiles
quadrat.test(pattern, nx = 4, ny = 4)
##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern
## X2 = 21.654, df = 15, p-value = 0.2344
## alternative hypothesis: two.sided
##
## Quadrats: 4 by 4 grid of tiles
Se trata de una estimación no paramétrica de la función de intensidad.
d <- density(pattern, edge = TRUE, kernel="gaussian")
plot(d, main = "Gaussian kernel density stimation GFP + RFP")
Este modelo asume que el proceso es estacionario (proceso estocástico con una distribución de probabilidad relativamente constante).
rescaled.pattern <- rescale(pattern, 10)
fryplot(rescaled.pattern, main = "Fry GFP + RFP")
La función empírica K asume asume de nuevo un proceso estacionario.
plot(Kest(pattern), main="K function")
Cada función corresponde a diferentes correcciones de los ejes: * theo Kpois(r) theoretical Poisson K(r)
* border Kbord(r) border-corrected estimate of K(r)
* trans Ktrans(r) translation-corrected estimate of K(r) * iso Kiso(r) isotropic-corrected estimate of K(r)
Transofrmando la función K a un alinea recta mediante la función L de Besag obtenemos la siguiente gráfica que nos permite apreciar mejor las desviaciones. La función aplica las mismas correcciones que aplicaba antes.
plot(Lest(pattern), main = "Besag’s L-function")
La función de correlacción por pares (pair correlation function) considera en el numerador la probabilidad de observar un par de puntos separados por una distancia r, dividida por esa misma probabilidad para un proceso estocástico. En este caso se puede ver que hay un claro patrón de clustering dentro de una pequeña distancia. Esto puede ser indicativo de posibles agrupaciones entre un número reducido de bacterias.
plot(pcf(pattern), main = NULL)
Se han seleccionado 99 simulaciones, de esta forma, la hipótesis nula es rechazada con un valor de alfa = 0.01.
a <- capture.output(plot(envelope(pattern, Kest, nsim = 99), main = NULL, legend = FALSE))
a <- capture.output(plot(envelope(pattern, Lest, nsim = 99), main = NULL, legend = FALSE))
La información referente al espaciado entre los puntos proporciona información complementaria al estudio de la correlación realizado anteriormente. La función empty-space (función F) asumiendo de nuevo un estado estacionario corresponde a la función de distribución empírica de las distancias observaas respecto al espacio vacio de m localizaciones. En este caso tenemos una distancia acumulativa. De nuevo la librería spatstats aplica diferentes correcciones.
plot(Fest(pattern), main = "F-function (point-to-event)")
a <- capture.output(plot(envelope(pattern, Fest, nsim = 99), main = NULL, legend = FALSE))
La función G (Nearest neighbour distances o event-to-event) de nuevo asumiendo un estado estacionario computa la probabildiad de que un punto tenga una distancia menor, mayor o igual a un proceso estacionario.
plot(Gest(pattern), main = "G-function (event-to-event)")
a <- capture.output(plot(envelope(pattern, Gest, nsim = 99), main = NULL, legend = FALSE))
La función J es una combinación de las funciones F y G.
plot(Jest(pattern), main = "J-function")
En este caso suponiendo que los datos no tengan un patrón homogéneo de intensidad.
plot(Linhom(pattern), main = "L-function")
plot(Jinhom(pattern), main = "J-function")
plot(Ginhom(pattern), main = "G-function")
plot(Finhom(pattern), main = "F-function")
a <- capture.output(plot(envelope(pattern, Linhom, nsim = 99), main = "L-function", legend = FALSE))
a <- capture.output(plot(envelope(pattern, Jinhom, nsim = 99), main = "J-function", legend = FALSE))
a <- capture.output(plot(envelope(pattern, Ginhom, nsim = 99), main = "G-function", legend = FALSE))
a <- capture.output(plot(envelope(pattern, Finhom, nsim = 99), main = "F-function", legend = FALSE))
La metodología seguida ha sido la misma pero aplicada de forma independiente a cada una de las bacterias.
green.data <- data.frame(x = image.df[image.df$value > 0.5 & image.df$cc == 2 & image.df$quality != 0, ]$x,
y = image.df[image.df$value > 0.5 & image.df$cc == 2 & image.df$quality != 0, ]$y,
strain = image.df[image.df$value > 0.5 & image.df$cc == 2 & image.df$quality != 0, ]$cc)
red.data <- data.frame(x = image.df[image.df$value > 0.5 & image.df$cc == 1 & image.df$quality != 0, ]$x,
y = image.df[image.df$value > 0.5 & image.df$cc == 1 & image.df$quality != 0, ]$y,
strain = image.df[image.df$value > 0.5 & image.df$cc == 1 & image.df$quality != 0, ]$cc)
p.green <- ggplot(green.data, aes(x=x, y=y)) +
geom_point(aes(color = as.factor(strain))) + theme_minimal() + theme(legend.title = element_blank()) +
scale_color_manual(values = c("#04FF00"))
p.red <- ggplot(red.data, aes(x=x, y=y)) +
geom_point(aes(color = as.factor(strain))) + theme_minimal() + theme(legend.title = element_blank()) +
scale_color_manual(values = c("#FF0000"))
plot_grid(p, p.green, p.red, nrow = 1, labels = c("A", "B", "C"))
pattern.green <- ppp(green.data$x, green.data$y, c(0, 112), c(0, 112))
pattern.green <- rescale(pattern.green, 10)
pattern.red <- ppp(red.data$x, red.data$y, c(0, 112), c(0, 112))
pattern.red <- rescale(pattern.red, 10)
intensity(pattern.green)
## [1] 5.030293
sqrt(intensity(pattern.green)/area(Window(pattern.green)))
## [1] 0.2002528
intensity(pattern.red)
## [1] 0.6776148
sqrt(intensity(pattern.red)/area(Window(pattern.red)))
## [1] 0.07349764
Q.green <- quadratcount(pattern.green, nx=6, ny=6)
Q.red <- quadratcount(pattern.red, nx=6, ny=6)
quadrat.test(pattern.green, nx = 6, ny = 6)
##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern.green
## X2 = 99.441, df = 35, p-value = 8.652e-08
## alternative hypothesis: two.sided
##
## Quadrats: 6 by 6 grid of tiles
quadrat.test(pattern.red, nx = 4, ny = 4)
##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern.red
## X2 = 38.671, df = 15, p-value = 0.001435
## alternative hypothesis: two.sided
##
## Quadrats: 4 by 4 grid of tiles
plot(intensity(Q.green, image = TRUE), main = "Intensity per quadrant in GFP")
plot(intensity(Q.red, image = TRUE), main = "Intensity per quadrant in RFP")
plot(density(pattern.green, edge = TRUE, kernel="gaussian"),
main = "Gaussian kernel density stimation in GFP")
plot(density(pattern.red, edge = TRUE, kernel="gaussian"),
main = "Gaussian kernel density stimation in RFP")
fryplot(pattern.green, main = "Fry GFP")
fryplot(pattern.red, main = "Fry RFP")
plot(Kest(pattern.green), main="K function GFP")
plot(Lest(pattern.green), main="L function GFP")
plot(Kest(pattern.red), main="K function RFP")
plot(Lest(pattern.red), main="L function RFP")
plot(pcf(pattern.green), main = "GFP")
plot(pcf(pattern.red), main = "RFP")
a <- capture.output(plot(envelope(pattern.green, Kest, nsim = 99), main = "K-function GFP", legend = FALSE))
a <- capture.output(plot(envelope(pattern.green, Lest, nsim = 99), main = "L-function GFP", legend = FALSE))
a <- capture.output(plot(envelope(pattern.red, Kest, nsim = 99), main = "K-function RFP", legend = FALSE))
a <- capture.output(plot(envelope(pattern.red, Lest, nsim = 99), main = "L-function RFP", legend = FALSE))
plot(Fest(pattern.green), main = "F-function GFP")
plot(Gest(pattern.green), main = "G-function GFP")
plot(Jest(pattern.green), main = "J-function GFP")
plot(Fest(pattern.red), main = "F-function RFP")
plot(Gest(pattern.red), main = "G-function RFP")
plot(Jest(pattern.red), main = "J-function RFP")
a <- capture.output(plot(envelope(pattern.green, Fest, nsim = 99), main = "F-function GFP"))
a <- capture.output(plot(envelope(pattern.green, Gest, nsim = 99), main = "G-function GFP"))
a <- capture.output(plot(envelope(pattern.green, Jest, nsim = 99), main = "J-function GFP"))
a <- capture.output(plot(envelope(pattern.red, Fest, nsim = 99), main = "F-function RFP"))
a <- capture.output(plot(envelope(pattern.red, Gest, nsim = 99), main = "G-function RFP"))
a <- capture.output(plot(envelope(pattern.red, Jest, nsim = 99), main = "J-function RFP"))